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ABSTRACT 



Aims. We extend our framework for 3D radiative transfer calculations with a non-local operator splitting methods along (full) charac- 
teristics to spherical and cylindrical coordinate systems. These coordinate systems are better suited to a number of physical problems 
than Cartesian coordinates. 

Methods. The scattering problem for line transfer is solved via means of an operator splitting (OS) technique. The formal solution is 
based on a full characteristics method. The approximate A operator is constructed considering nearest neighbors exactly. The code is 
parallelized over both wavelength and solid angle using the MPl library. 

Results. We present the results of several test cases with different values of the thermalization parameter for the different coordinate 
systems. The results are directly compared to ID plane parallel tests. The 3D results agree very well with the well-tested ID calcula- 
tions. 

Conclusions. Advances in modern computers will make realistic 3D radiative transfer calculations possible in the near future. 
Key words. Radiative transfer - Scattering 



1. Introduction 



In iHauschildt & Baroiil ilOOd 



hereafter: Paper 
hereafter: Paper II) 



I), 
and 



[Baron & Hauschildt' (TOOTI 
(iHauschil dt & Baron 2008- hereafter: Paper III) we described a 
framework for the solution of the radiative transfer equation for 
scattering continua and lines in 3D (when we say 3D we mean 
three spatial dimensions, plus three momentum dimensions) for 
the time independent, static case using full characteristics based 
on voxels in Cartesian coordinates. Cartesian coordinates are 
by no means an optimal choice for many 3D physical systems, 
for example supernova atmospheres are better described in 
3D spherical and accretion disks are better described in 3D 
cylindrical coordinates. These coordinate systems make it a 
little harder and computationally more expensive to track the 
characteristics through the curvilinear voxel grid, however, 
this is more than overcompensated by the higher efficiency of 
voxel usage, for example consider putting a roughly spherical 
structure where density and composition vary strongly with 
radius into a Cartesian coordinate system cube. 

We describe our method, its rate of convergence, and present 
comparisons to our well-tested 1-D calculations. 



are discussed in detail in Papers I and II and will not be repeated 
here. 

We use a "full characteristics" approach that completely 
tracks a set of characteristics of the radiative transfer equation 
from the inner or outer boundary through the computational do- 
main to their exit voxels and takes care that each voxel is hit 
by at least one characteristic per solid angle. One characteristic 
is started at each boundary voxel and then tracked through the 
voxel grid until it leaves the other boundary. The direction of a 
bundle of characteristics is determined by a set of solid angles 
(0, (f) which correspond to a normalized momentum space vec- 
tor (pj,, py, P-). 

In contrast to the case of Cartesian coordinates the volumes 
of the voxels can be vastly different in spherical and Cartesian 
coordinate systems. Therefore, tracking a characteristic across 
the voxel grid is done using an adaptive stepsize control along 
the characteristic so that a characteristic steps from one voxel to 
its neighbor in the direction of the characteristic. In addition, we 
make sure that each voxel is hit by at least one characteristic for 
each solid angle, which is important for the inner (small radii) 
voxels. 

The code is parallelized as described in Papers I and II. 



2. Method 

In the following discussion we use notation of Papers I - III. The 
basic framework and the methods used for the formal solution 
and the solution of the scattering problem via operator splitting 



3. spherical coordinates 

3.1. Testing environment 

We use the framework discussed in Paper I and II as the baseline 
for the line transfer problems discussed in this paper Our basic 
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setup is similar to that discussed in Paper II. We use a sphere 
with a grey continuum opacity parameterized by a power law in 
the continuum optical depth Tstd- The basic model parameters are 



1. Inner radius = 10 cm, outer radius rout = 2 x 10 cm 
(large test case) or rout = 1-1 x 10'" cm (small test case). 

2. Minimum optical depth in the continuum t™™ = lO"'' and 
maximum optical depth in the continuum t™^" = 10*^. 

3. Constant temperature with Teff - 10"^ K, therefore, the 
Planck function B is also constant. 

4. Outer boundary condition I^^ s and diffusion inner bound- 
ary condition for all wavelengths. 

5. Continuum extinction ;ife = C/r^, with the constant C fixed 
by the radius and optical depth grids, for the continuum tests, 
Xc = C for the line transfer tests 

6. Parameterized coherent & isotropic continuum scattering by 
defining 

Xc = ecKc + (1 - ec)o-c 

with < e,. < I. Kc and cTc are the continuum absorption and 
scattering coefficients. 

This problem is used because the results can be directly com- 
pared with the results obtained with a ID spherical radiation 
transport code (IHauschildti 1 1 992b to assess the accuracy of the 
method. The sphere is put at the center of the 3D grid. The radial 
grid in the 3D coordinates is logarithmic in Tstd and has, there- 
fore, highly variable spacing in r. In 3D spherical coordinates 
(r, 6c, 4>c) (note: the subscript c indicates a coordinate in order to 
distinguish with the solid angle Q = {6, (p)) we use regular grid 
in {0c, (pc) with resolutions given below. In the tests with cylin- 
drical coordinates (p, (pc, z) we use a regular grid in z and (pc and 
a variable grid in p. The solid angle space was discretized in 
(0, (p) with «fl - if not stated otherwise. In the following we 
discuss the results of various tests. In all tests we use the full 
characteristics method for the 3D RT solution. 

The line of the simple 2-level model atom is parameterized 
by the ratio of the profile averaged line opacity xi to the contin- 
uum opacity Xc and the line thermalization parameter e/. For the 
test cases presented below, we have used;^'^ = const, ec - I, and 
a grey temperature structure. We set xi/Xc - 10^ ^ simulate a 
strong line, with varying e/ (see below). With this setup, the op- 
tical depths as seen in the line range from 10"~ to 10*. We use 32 
wavelength points to model the full line profile, including wave- 
lengths outside the line for the continuum. We did not require 
the line to thermalize at the center of the test configurations, this 
is a typical situation for a full 3D configuration as the location 
(or even existence) of the thermalization depths becomes more 
ambiguous than in the ID case. 

In the following we discuss the results of various tests. 
Unless otherwise stated, the tests were run on parallel computers 
using 4 to 5 12 CPUs (depending on the test). In spherical coordi- 
nates we use a voxel grid with (rir, ng^, n^J - (65, 33, 65) voxels, 
this corresponds to an angular resolution of about 5.5 degrees in 
0c and (pc- Projected on the surface of the Earth, this means that 
Europe is covered by about 43 voxels. The spherical test cases 
use incoming intensities of zero at the outside (Rm-dx) and I - B 
at the inside {{Rmm)- 



3.2. Results 

3.2.1. Continuum tests 

In Fig. [T] we show the results for 3D continuum transfer with 
3D spherical coordinates for various values of the thermalization 
parameter e^. . The results are compared to ID results (full line). 
In all cases, the agreement is excellent, even for the case with 
ec = 10"*. Compared to paper I the agreement is better, this is 
caused by the better discretization of a spherical configuration in 
spherical coordinates, rather than fitting a sphere into a Cartesian 
grid. It illustrates the importance of adapted coordinate systems 
to obtain the optimal mix between computational efficiency and 
numerical accuracy. 

Good coverage of the solid angle is important for 3D transfer. 
The ID code automatically adjusts the angle coverage related to 
the geometry of the conffguration, however, in the 3D code the 
solid angles have to be prescribed. In Fig.|2]we demonstrate the 
importance of solid angle sampling for a 'small' test case with 
fc = 10"^. The plots show the ratio J/B as function of radial 
optical depth for various solid angle samplings and for all {(pc, dc) 
compared to the ID results. It is immediately clear that for the 
'small' test case a solid angle sampling of 64^ points or more 
is reasonable, but that 16^ or less points result in a considerable 
scatter of the results. This result does not depend strongly on the 
radial extension of the conffguration, as shown in Fig. [3] 

The convergence rate of the nearest neighbor A* (3^ neighbor 
voxels are considered in the A*, see Paper I for details) with Ng 
acceleration is virtually identical to that of the ID tri-diagonal 
A* case, cf. Fig. |4] The importance of a non-local A* for good 
convergence in the 3D case can be seen by comparing the re- 
sults obtained with a local (diagonal) A*. Its convergence with 
Ng acceleration is actually worse than that of the non-local A* 
without Ng acceleration. The excellent convergence of the non- 
local A* can be used to reduce the computing time by using a 
less stringent convergence criterion that still gives results accu- 
rate enough for a particular application. One has to be careful to 
be not too aggressive in relaxing the convergence criterion. This 
is highlighted in Fig.|5]where the results for the large test case are 
shown for different combinations of A* with Ng acceleration and 
difference convergence limits. It is obvious that the combination 
of Ng acceleration and a relatively large convergence limit with 
a local (diagonal) A* is not acceptable, whereas the same con- 
vergence limit (10"^ maximum relative change of /) with a non- 
local (nearest neighbor) A* gives usable results even without Ng 
acceleration. Setting a convergence limit of 10""* gives reason- 
able results in all test cases and can reduce computing time by 
about 50%. 

With 128 MPI processes on the SGI ICE system of the 
Hochstleistungs Rechenzentrum Nord (HLRN) the continuum 
test with e - 10"'* and = - 65 and iig - 33 spatial points 
and 256^ solid angle points requires 16 iterations and 1543 s 
wallclock time, where a formal solution takes about 80s, the 
construction of the non-local A* about 140 s and about 1.6 s are 
used for the solution of the operator splitting step. For the same 
case with 64- solid angles the total wallclock time is 230 s. The 
memory requirements are about 45MB/process. 
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3.2.2. Line tests 

In Fig.|6]we show the resuhs of a comparison of the ratio of the 
line mean intensities / and the Planck function B for the 'large' 
test case and various values of e/. The ID comparison models 
assume a diffusive inner boundary condition, whereas the 3D 
models use I - B as the boundary condition for the inner bound- 
ary. Therefore, the models will differ close to the inner bound- 
ary. Away from the inner boundary the results of the ID and 3D 
calculations agree very well with each other. The number of it- 
erations required to reach a prescribed accuracy is similar to the 
continuum tests described above. This shows that the quaUty of 
the 3D calculations match the ID results, which is very impor- 
tant for future applications. 

In Figs. |7] and [8] we show the flux vectors of the outermost 
voxels (outermost radius, all {6c, (pc) points) for the line transfer 
case with e/ = 1 and for different numbers of solid angle points. 
In the test case, the radial component of the flux vector Fr should 
be the only non-zero component as in spherical symmetry both 
Fg and F^ are zero. The figures show that this is indeed the case 
if a large enough number of solid angle points is used in the cal- 
culation, if the number is too small, Fg and F^ can have errors of 
about 10% and F,. scatters considerably around the ID compar- 
ison. For 512^ solid angle points, the scatter of F,. is too small 
to be seen in the plot and Fg and F^ are smaller than 1 % of F,- 
over the whole surface of the spherical grid. This highlights the 
necessity of having enough solid angle points to obtain good ac- 
curacy. Fig.|9]shows Fr compared to the ID results for a number 
of values of ei. In all cases the results agree very well with each 
other, showing the good numerical accuracy of the 3D code. The 
results of these tests are valid only for the test cases, therefore, 
a similar test should be performed for each "real" 3D configura- 
tion individually in order to make sure that enough solid angles 
are used for the accuracy requirements of each individual prob- 
lem. 

With 512 MPI processes on the SGI ICE system of the 
HLRN the line test with e/ - 10 and - «0 = 65 and 
rig = 33 spatial points, 32 wavelength points and 256^ solid an- 
gle points requires 30 iterations and 20287 s wallclock time and 
for the same case with 64- solid angles the total wallclock time 
is 1400 s. In both cases we used 32 "clusters" with 16 MPI pro- 
cesses each so that for each wavelength 16 MPI processes work 
on a formal solution (see paper II for details). 

4. Cylindrical coordinates 

The setup for the test of the code in cylindrical coordinates is es- 
sentially identical to the tests of the spherical coordinate system 
version. However, describing a spherically symmetric objects 
in cylindrical coordinates is non-optimal and we expect results 
comparable to those presented in papers I and II. The cylindrical 
coordinate system version of the code allows for different maps 
in p and z for flexibility describing objects with nearly cylindri- 
cal symmetry. 

The differences in the basic setup compared to the spherical 
coordinates tests are 

1. inner radius rc - 10'° cm, outer radius rout = 11 x lO'^cm. 

2. linear radial grid. 



3. Continuum extinction ;\fc = C/r^, with the constant C fixed 
by the radius and optical depth grids, for the continuum tests, 
Xc - C for the line transfer tests. 

All other test parameters are the same as above. 
4.1. Results 

For the cylindrical coordinate tests we show just the results of 
a standard test discussed in Paper I and above for the spheri- 
cal case with e - 10"^. The results for the continuum test case 
are presented in Fig. [TOl The results agree reasonably well with 
the ID test case, comparable to the results described in Paper 
I. The reason for this is that describing a sphere in cylindrical 
coordinates is not an ideal use of this coordinate system, there- 
fore much higher grid resolution is required to reach reasonable 
accuracy (the results for grids with «p = iig = «; = 64 are sub- 
stantially worse). As for the test cases in spherical coordinates 
the resolution in solid angles is very important for accurate re- 
sults. The convergence behavior is the same as in the case of 
spherical coordinates. 

5. Conclusions 

We have discussed the extension of our 3D framework to 3D 
spherical and cylindrical coordinate systems. These coordinates 
are better suited than Cartesian coordinates for problems that are 
approximately spherical (e.g., supernova envelopes) or cylindri- 
cal (e.g., protoplanetary disks). The convergence properties and 
the accuracy of the solutions are comparable to the ID solution 
for the test cases discussed here. These additional coordinate 
systems allow for more general use with higher internal accu- 
racy while using less memory than the purely Cartesian coordi- 
nate system. 

The next step in the development of the 3D framework is to 
add the treatment of velocity fields in the observer's frame and 
the co-moving frame. The latter will be important for applica- 
tions with relativistic velocity fields, e.g., supernovae, while the 
former is interesting for modeling the spectra of convective flow 
models. 
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Fig. 1. The ratio of the mean intensity J and the Planck function 
B as function of the radial optical depth Tstd for the 'large' test 
case in spherical coordinates and different values e^. . The optical 
depth is measured along the radius r. 
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Fig. 2. Effect of solid angle resolution on the results of the 3D 
calculations. The plots show the results for the 'small' contin- 
uum test case in spherical coordinates with ec - 10 and for 
various solid angle resolutions with vg - from (256, 256) to 
(8, 8). Full line; results of the ID calculation, dots; results from 
the 3D calculations. 



Fig. 3. Effect of solid angle resolution on the results of the 3D 
calculations in spherical coordinates. The plots show the results 
for the 'large' continuum test case with e^. - 10 and for various 
solid angle resolutions with vg - from (256, 256) to (8, 8). 
Full line; results of the ID calculation, dots; results from the 3D 
calculations. 
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Fig. 4. Convergence behavior for the 'large' continuum test case 
in spherical coordinates with ec - 10""^ for various combinations 
of local/non-local A* and Ng acceleration. 'NN' indicates the 
non-local nearest neighbor A* operator, 'diagonal' the local A*. 
For comparison the results for the ID spherical code and the 3D 
lambda iterations are also plotted. 
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Fig. 5. Effects of relaxing the convergence criterions on the re- 
sults of the 3D calculations in spherical coordinates. The plots 
show the results for the Targe' continuum test case with = 
10""^ and for various convergence limits for nearest neighbor 
(NN) and diagonal (Diag) A* operators with and without Ng 
acceleration. The numbers are the maximum allowed relative 
change over all voxels used to halt the iterations. Full line: results 
of the ID calculation, dots: results from the 3D calculations. 
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Fig. 7. Flux vectors in the outermost voxels for the line trans- 
fer test in spherical coordinates with e/ = 1 and 64^ solid angle 
points. The top panel shows the radial component Fr of the voxel 
flux vector compared to the results for the ID comparison calcu- 
lation. The middle and bottom panels show the ratio of the polar 
Fe and azimuthal F^ flux components to Fr, respectively. 



Fig. 6. The ratio of the line mean intensity / and the Planck 
function B as function of the radial optical depth Tstd for the 
Targe' test case in spherical coordinates for different values of 
€[. The optical depth is measured along the radius r. 
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Fig. 10. Results for an example of a grey radiative transfer cal- 
culation for a sphere in cylindrical 3D coordinates for e - 10"^ 
with Hp - Tig = «- - 64 and 64^ solid angle points. The results of 
the ID calculations are shown for comparison with * symbols. 
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Fig. 8. Flux vectors in the outermost voxels for the line trans- 
fer test in spherical coordinates with ei - I and 512^ solid angle 
points. The top panel shows the radial component Fr of the voxel 
flux vector compared to the results for the ID comparison calcu- 
lation. The middle and bottom panels show the ratio of the polar 
Fg and azimuthal F^ flux components to Fr, respectively. 
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Fig. 9. Radial components Fr of the flux vectors in the outermost 
voxels for the line transfer tests in spherical coordinates with 
ei = 1, 10"^, 10"^, 10"** and 256^ solid angle points. The results 
of the ID calculations are shown for comparison. 



